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ABSTRACT 

This paper proposes a new semi-analytic modelling of galaxy properties in the 
IR/submm wavelength range, which is explicitly set in a cosmological framework. We 
start from a description of the non-dissipative and dissipative collapses of primordial 
perturbations, and add star formation, stellar evolution and feedback, as well as the 
absorption of starlight by dust and its re-emission in the IR and submm. This type of 
approach has had some success in reproducing the optical properties of galaxies. We 
hereafter propose a simple extension to the IR/submm range. The growth of struc- 
tures is followed according to the standard Cold Dark Matter model. We assume that 
star formation proceeds either in a "quiescent" mode, e.g. as in disks, or in a "burst" 
mode with ten times shorter time scales. In order to reproduce the current data on 
the evolution of the comoving cosmic SFR and gas densities, we need to introduce a 
mass fraction involved in the "burst" mode strongly increasing with redshift, probably 
reflecting the increase of interaction and merging activity. We estimate the IR/submm 
luminosities of these "mild starburst" and "luminous UV/IR galaxies" , and we explore 
how much star formation could be hidden in heavily-extinguished, "ultraluminous IR 
galaxies" by designing a family of evolutionary scenarios which are consistent with 
the current status of the "cosmic constraints" , as well as with the IRAS 60 /im lu- 
minosity function and faint counts, but with different high-z IR luminosity densities. 
However, these scenarios generate a Cosmic Infrared Background whose spectrum falls 
within the ±lcr range of the isotropic IR component detected by Puget et al. (1996) 
and revisited by Guiderdoni et al. (1997). We give predictions for the faint galaxy 
counts and redshift distributions at IR and submm wavelengths. The submm range is 
very sensitive to the details of the evolutionary scenarios. As a result, the on-going 
and forthcoming observations with ISO and SCUBA (and later with SIRTF, SOFIA, 
FIRST and PLANCK) will put strong constraints on the evolution of galaxies at z ~ 1 
and beyond. 
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1 INTRODUCTION 

This paper describes a new modelling of galaxy evolution in 
the infrared and submm wavelength ranges, which are now 
open to high-redshift exploration by the on-going observa- 
tions with ISO and SCUBA, and forthcoming facilities and 
experiments such as SIRTF, SOFIA, FIRST and PLANCK. 

As a matter of fact, our knowledge of the early epochs 
of galaxies has recently increased thanks to the richness 
and precision of the observational evidence obtained by 
UV/visible/NIR surveys of high-redshift objects (Lilly et al. 
1995; Ellis et at 1996; Cowie et al 1996; Steidel et al. 1996; 
Williams et al. 1996). The pattern of galaxy evolution which 



emerges from these data can be summarized as follows: i) 
Faint galaxy counts show the presence of a large number 
of blue objects, well in excess of no-evolution expectations 
(Williams et al. 1996). ii) These blue objects are "sub-L*" 
galaxies undergoing strong bursts of star formation (Lilly et 
al. 1996; Ellis et al. 1996). iii) The fraction of these blue ob- 
jects with unclassified/peculiar morphologies showing signs 
of tidal interaction and merging (Abraham et al. 1996) in- 
creases from local samples to the HST high-resolution ob- 
servations of the Medium Deep Survey (Griffiths et al. 1994) 
and Hubble Deep Field (Williams et al. 1996). iv) The global 
star formation rate (hereafter SFR) of the universe declined 
by a factor of about ten since redshift z ~ 1 (Lilly et al. 
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1996; Madau et al. 1996; Sawicki et al. 1997; Connolly et 
al. 1997). v) Finally, this high SFR seems to be correlated 
to the decrease of the cold-gas comoving density associated 
with damped Lyman-a systems between z — 2 and z = 
(Storrie-Lombardi et al. 1996). These results nicely match a 
picture in which star formation in bursts triggered by inter- 
action/merging consumes the gas content of galaxies as time 
goes on. It is common wisdom that such a qualitative sce- 
nario is expected within the paradigm of hierarchical growth 
of structures. The implementation of hierarchical galaxy for- 
mation in semi-analytic models quantitatively substantiates 
this view and may further suggest that we have seen the 
bulk of star formation and understood the broad features of 
galaxy formation (Baugh et al. 1997). 

However, it should be emphasized that this seemingly 
consistent view is entirely based on visible/NIR observa- 
tions which only probe the rest-frame UV/visible light of 
high-z objects. The total amount of energy released by 
star formation should be estimated by summing up the 
UV/visible/NIR light of stellar populations directly escap- 
ing from galaxies, and the part which has been absorbed by 
dust and re-emitted in the IR/submm wavelength range. 
The corrections needed to account for optical extinction by 
dust are rather uncertain and could easily induce an upward 
revision of the high-redshift SFR deduced from rest-frame 
UV/visible observations by factors of a few. Moreover, a sig- 
nificant fraction of star formation might be completely hid- 
den in heavily-extinguished galaxies which are missed by 
the above-mentionned surveys. Even for normal objects like 
spiral galaxies, the fraction of energy emitted in the IR can 
amount to 30 % of the optical, and this can increase to more 
than 95 % for ultraluminous starbursts. As a matter of fact, 
the IR/submm range might be more adequate than the op- 
tical for observing high-redshift, starburst galaxies provided 
they are dusty. The rest-frame IR/submm spectrum of dust 
heated by a young stellar population peaks between 50 to 
100 /im, pointing out the submm range as particularly rele- 
vant for the search of primeval, high-redshift galaxies. 

At this time, we know very little about the "optically- 
dark" side of galaxy evolution. The IRAS satellite has given 
the first complete survey of FIR galaxy properties, in four 
bands between 12 and 100 fj,m. Many studies have empha- 
sized the wide variety of IR luminosities, from "normal spi- 
rals" to "mild starbursts", and finally to the "luminous 
IR galaxies" (hereafter LIGs), mostly interacting systems, 
and the spectacular "ultraluminous IR galaxies" (hereafter 
ULIGs), which are mergers (Sanders and Mirabel 1996; 
Clements et al. 1996b). Unfortunately, while we have learned 
from IRAS that about one third of the bolometric luminos- 
ity of the local universe is released in the IR/submm (Soifer 
and Neugebauer 1991), we know very little about galaxy 
evolution in this wavelength range. Faint galaxy counts and 
redshift surveys down to flux densities S v ~ 60 mjy (at 
60 (im) do not probe deeper than z ~ 0.2 (Ashby et al. 
1996; Clements et al. 1996a). These surveys seem to show a 
strong luminosity and/or density evolution of IRAS sources, 
but it is difficult to extrapolate this trend to higher redshifts 
on a firm ground. In spite of its observational limits, IRAS 
has also revealed the existence of IRAS 10214+4724, a very 
peculiar, "hyper luminous" galaxy at z = 2.286 (Rowan- 
Robinson et al. 1991a), though this object is likely affected 
by lensing (Eisenhardt et al. 1996). It is expected that the 
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Figure 1. Observer-frame model spectra of a L/^ = 10 Lq 
galaxy at increasing redshifts (from top to bottom), for a cos- 
mology with h = 0.5 and f2o = 1. The details of the modelling 
are explained in Sect. 3.2. The reader is invited to note that the 
apparent flux in the submm range is almost insensitive to red- 
shift, because the shift of the 100 (im bump counterbalances the 
distance dimming. 



ISO satellite will complete and detail this picture, in a 
broader wavelength range from a few ftm to 200 /im. Other 
projects, such as SIRTF, will give access to better sensitivity 
and imaging capabilities. 

The properties of galaxies in the submm range are sen- 
sitive to the spectral characteristics of dust, especially its 
cmissivity at large wavelengths which is not constrained 
by IRAS observations alone. With respect to the relative 
wealth of data in the FIR, the submm emission of galaxies 
is poorly known. The observational literature gathers the 
submm fluxes of only a few tens galaxies which have been 
measured from ground-based or aircraft-borne instruments. 
These observations are difficult and some of the estimates 
of the amount of energy released in this range happen to 
be strongly discrepant (e.g. Chini et al. 1986; Stark et al. 
1989; Eales et al. 1989; Chini and Kriigel 1993). However, 
the observational situation will soon evolve with the start of 
SCUBA operations (and later, with SOFIA and especially 
FIRST). About twenty counterparts of radiogalaxies and 
quasars have already been observed at submm/mm wave- 
lengths (see e.g. Hughes et al. 1997; and references therein). 
Although their evolutionary status is not fully understood, 
these objects give a first flavour of the future deep surveys 
which will be achieved by forthcoming instruments. 

Theoretical modelling is needed to optimise observa- 
tional strategies and also to help assess the point source 
contribution of the forthcoming CMB satellites MAP, and 
especially PLANCK. These experiments aim at mapping the 
anisotropies of the submm/mm sky on scales well below the 
degree. On these scales, the separation of the various fore- 
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Figure 2. Number density and formation rate of collapsed haloes, 
and their relation to two modes of star formation. The left-hand 
panels show the number density of collapsed haloes (top) and 
the formation rate of collapsed haloes (bottom) computed for the 
SCDM model with h = 0.5, H = 1 and tr 8 = 0.67. The curves 
are plotted for redshifts z = 3.67, 1.99, 0.92, 0.23 and (for 
increasing number densities at the high-mass end). The right- 
hand top panel shows the evolution of the IR luminosity function 
at the same redshifts, in the so-called "quiescent" mode of star 
formation with /3 = 100, a = 5 and V^oi = 130 km s _1 . The 
increase of the number of galaxies at all masses, as time goes on, 
is similar to the increase of the number of haloes and reflects the 
accumulation of galaxies. The right-hand bottom panel shows the 
evolution of the IR luminosity function at the same redshifts, in 
the so— called "burst" mode of star formation with /3 = 10. The 
galaxies undergo strong starbursts which act as beacons at the 
epochs of their initial collapses. Consequently, the evolution of the 
luminosity function reflects the formation rate of new haloes, with 
the characteristic crossing of the luminosity functions between low 
and high masses. 



grounds and backgrounds which superimpose to the fluctu- 
ations of the CMB is more difficult than for the 7° field 
of view of the very successful COBE/~DMR experiment. In 
the case of PLANCK, preliminary studies have shown that 
several thousands of galaxies should be detected at several 
wavelengths between 400 /im and 1 cm. 

The epoch of galaxy formation can also be observed by 
its imprint on the background radiation which is produced 
by the line-of-sight accumulation of extragalactic sources. 
The search for the "Cosmic Optical Background" (hereafter 
COB) currently gives only upper limits. Nevertheless, the 
shallowing of the faint counts obtained in the HDF sug- 
gests that we are now close to convergence (Williams et al. 
1996). Thus the lower limit to the COB obtained by sum- 
ming up the contributions of faint galaxies is likely close to 
the real value. At longer wavelengths, the DIRBE instru- 
ment on COBE has given upper limits on the FIR back- 
ground between 2 and 300 /im (Hauser 1995), while Puget 



Figure 3. Distribution of the characteristic SFR time scales for 
the "quiescent" mode of star formation (/3 = 100, Vhot = 130 km 
s _1 and a = 5). Only galaxies with V c > are retained. The 

histogram shows the data for a sample of 63 bright disk galaxies 
analysed by Kennicutt et al. (1994). 

et al. (1996) have discovered an isotropic component in the 
COBE/FIRAS residuals between 200 /im and 2 mm, which 
could be the long-sought "Cosmic Infrared Background" 
(hereafter CIB). The presence of this component is con- 
firmed by a new analysis restricted to the cleanest regions 
of the sky, where the foreground Galactic components are 
negligible (Guiderdoni et al. 1997). Such a detection seems 
to yield the first "post-Z/MS" constraint on the high-z evo- 
lution of galaxies in the IR/submm range, before the era of 
ISO results. Its level comparable to the above-mentionned 
estimate of the COB suggests that a significant fraction of 
the energy of young stars is absorbed by dust and released 
in the IR/submm. 

The models which have been proposed to predict the 
faint counts and background radiation in the IR/submm 
can be classified as "backward evolution" and "forward evo- 
lution", following the good review by Lonsdale (1996). In 
the first class of models, the IR/submm luminosity function 
undergoes luminosity and/or number evolution which are 
simply parameterized as power laws of (1 + z) (e.g. Weed- 
man 1990; Beichman and Helou 1991; Hacking and Soifer 
1991; Oliver et al. 1992; Treyer and Silk 1993; Blain and 
Longair 1993; Pearson and Rowan-Robinson 1996). These 
power laws are generally derived from fits of IRAS galaxy 
counts (which do not probe deeper than z ~ 0.2). Then 
they are extrapolated up to redshifts of a few units. Un- 
fortunately, various analyses of IRAS deep counts yield dis- 
crepant results at Seo < 300 mjy, and the issue of the evo- 
lution is very controverted (see e.g. Bertin et al. 1997, for a 
new analysis and discussion). 

In the second class of models, the photometric evolu- 
tion in the IR/submm is computed by implementing some 
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of the involved physical processes. For instance: chemical 
evolution which rules the amount of dust responsible for 
the IR/submm emission is modelled in Wang (1991a,b), 
and Eales and Edmunds (1996a,b). In addition to chemical 
evolution, the photometric evolution of stellar populations 
which heat the dust is modelled by Franceschini et al. (1991, 
1994) and Fall, Chariot and Pei (1996). While the previous 
models assume a simple relation between the dust content 
and the heavy-element abundance of the gas, Dwek and 
Varosi (1996) try to explicitly model the processes of dust 
formation and destruction. However, both classes of mod- 
els assume that all galaxies form at the same redshift Zf or 
and that there is no number evolution. But the paradigm of 
the hierarchical growth of structures implies that there is no 
clear-cut redshift z/ or since galaxy formation is a continuous 
process. Only Blain and Longair (1993a, b) proposed a for- 
malism to compute the redshift range of galaxy formation, 
in addition to chemical evolution, 

A consistent approach to the early evolution of galaxies 
is particularly important for any attempt at predicting their 
submm properties. Fig. 1 shows model spectra of a luminous 
IR galaxy as it would be observed if placed at different red- 
shifts. There is a wavelength range, between ~ 600 /im and 
~ 4 mm, in which the distance effect is counterbalanced 
by the "negative k-correction" due to the huge rest-frame 
emission maximum at ~ 100 /im. In this range, the apparent 
flux of galaxies depends weakly on redshift to the point that, 
evolution aside, a galaxy might be easier to detect at z = 5 
than at z — 0.5 ! The observer-frame submm fluxes, faint 
galaxy counts and diffuse background of unresolved galaxies 
are consequently very sensitive to the early stages of galaxy 
evolution. Note that this particular wavelength range brack- 
ets the maximum of emission of the CMB. 

The dissipative and non-dissipative processes ruling 
galaxy formation in dark matter haloes have been studied by 
various authors (see especially White and Rees 1978; Scha- 
effer and Silk 1985; Evrard 1989; Cole 1991; Blanchard et al. 
1992). The modelling of these processes, complemented by 
star formation, stellar evolution and stellar feedback to the 
interstellar medium, has been achieved at various levels of 
complexity, in the so-called semi-analytic approach which 
has been successfully applied to the prediction of the sta- 
tistical properties of galaxies (White and Frenk 1991; Lacey 
and Silk 1991; Lacey et al. 1993, Kauffmann et al. 1993, 
1994; Cole et al. 1994; Heyl et al. 1995; Kauffmann 1995, 
1996; Baugh et al. 1996a, b, 1997). It turns out that, in spite 
of differences in the details of the models, these studies lead 
to conclusions in the UV, visible and (stellar) NIR which are 
remarkably similar. 

None of these models has been applied so far to the pre- 
diction of the properties of galaxies in the IR and submm 
ranges. This is the main aim of this paper which proposes 
a simple version of the semi-analytic approach and gives 
predictions of faint counts at wavelengths between 60 pm 
and 1.4 mm. This first study also intends to identify some 
of the difficulties arising in such a modelling. Sec. 2 quickly 
reviews the main physical processes which have to be in- 
troduced in order to describe the formation of galaxies: the 
non-dissipative collapse of perturbations, the dissipative col- 
lapse of the gas, star formation and evolution, and feedback. 
Sec. 3 addresses the peculiar issue of dust absorption and re- 
emission in the IR and submm. Sec. 4 extracts useful infor- 
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Figure 4. Evolution of various quantities for a grid of models 
with ti, = 0.33, 1, 2, 3, 5, 10 and 20 Gyr. Panel a: gas fraction 
(from bottom to top). Panel 6: gas mctallicity (from top to bot- 
tom). Panel c: B — V colour (from top to bottom, taking into ac- 
count face-on internal extinction). Panel d: luminosity re— emitted 
by dust in the IR/submm (from top to bottom at 1 Gyr). 

mation from the recent UV/ visible deep surveys and uses the 
new constraint arising from the CIB in order to generate a 
family of evolutionary scenarios. Sec. 5 gives the IR/submm 
counts predicted from these scenarios. Finally, Sec. 6 briefly 
discusses these results, as well as the shortcomings of such 
models, and concludes. 

In a previous study, we showed that the CIB can be 
disentangled with a family of evolutionary scenarios which 
predict steep submm counts, and that the observations with 
ISO (at 175 /xm) will soon constrain the IR/submm evolu- 
tion at z ~ 1 and beyond (Guiderdoni et al. 1997). Here we 
present the details of the modelling and apply it to other ob- 
servations, more specifically in SCUBA bands (through the 
narrow submm atmospheric windows). A preliminary ver- 
sion of this work was presented in Guiderdoni et al. (1996). 
A forthcoming paper (Hivon et al. 1997) will show simula- 
tions of the anisotropies of the diffuse submm background 
due to galaxies and will study the possibility of their de- 
tection, with current and forthcoming instruments. Other 
papers in this series will try to overcome some of the short- 
comings of this first study. 



2 A SCHEMATIC VIEW OF GALAXY 
FORMATION 

2.1 Non-dissipative collapse 

The formation and evolution of a galaxy in its dark matter 
halo can be briefly sketched as follows: the initial pertur- 
bation, which is gravitationally dominated by non-baryonic 
dark matter, grows and collapses. After the (non-dissipative) 
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dots: data Franceschini & Andreani 1995 




Figure 5. Distribution of mean column densities at z = 0, in 
atom cm~ 2 , for galaxies with SB e < 25 mag arcsec -2 , and the 
"quiescent" mode of star formation ({3 = 100, a = 5 and Vjj t = 
130 km s -1 ). Dotted line: log < N to t >; solid line: log < N H >. 



collapse, and subsequent violent relaxation, the halo virial- 
izes, through the formation of a mean potential well seen 
by all particles, which consequently share the same velocity 
distribution. As detailed in the Appendix, we use the classi- 
cal top-hat model for spherically-symmetric perturbations, 
and we compute the mass distribution of collapsed haloes 
from the peaks formalism (Bardeen et al. 1986), as in Lacey 
and Silk (1991) and Lacey et al. (1993). Hereafter, we shall 
consider the SCDM model with H = 50 km s -1 Mpc -1 , 
0,0 = 1, A = 0, and erg = 0.67 as an illustrative case. We 
take a baryonic fraction Sls ar = 0.05, consistent with primor- 
dial nucleosynthesis. The redshift evolution of the number 
density and formation rate of collapsed haloes is plotted in 
fig. 2. 



2.2 Cooling and dissipative collapse 

The baryonic gas cools in the potential wells of dark matter 
haloes, by a process identical to cooling flows observed at 
the centre of rich clusters. The cooling time at halo radius 



. , , 3 n to t{r)kTv „ 2 2 



2n!(r)A(7V) 



Qbar A(Ty) ' 



(1) 



with /i e = 1.14 for ionized primordial gas. The cooling curve 
A(T) takes into account various cooling processes. We do not 
consider the metallicity dependence of A(T). Neglecting this 
dependence leads to an overestimate of cooling times, which 
are already very short. The equation t coo i(r coo i) = t(z) de- 
fines a cooling radius r coo i as a function of redshift z. At 
this redshift, only gas inside r cooi (or rv if r coo i > rv) 
cools and is available for star formation. This cooling cri- 
terion introduces a high-mass cut-off in the mass distribu- 



Figure 6. Distribution of mean, face-on optical depth in the 
B— band at z = 0, for the "quiescent" mode of star formation 
(/3 = 100, a = 5 and V^ ot = 130 km s~ 1 ). Two ways of comput- 
ing rg are considered. Dashes: Mean optical depth inside the gas 
radius (r g = 1.6r2s). Solid line: Mean optical depth inside one 
third of the Holmberg radius (at 26.5 mag arcsec -2 ). The dotted 
line shows the observational distribution in a sample of normal 
spirals and luminous IR galaxies from Franceschini and Andreani 
(1995). The total dust quantity derived from 1.3 mm observa- 
tions is distributed within one third of the Holmberg radius and 
an observational estimate of tb is computed. The observational 
distribution peaks at the value < rg >~ 1 and is brackctted by 
the two predicted histograms. 



tion of cold baryonic cores. At the low-mass end, the cool- 
ing is so efficient that almost all the gas can cool, lead- 
ing to a slope of the mass distribution of baryonic cores 
n(Mbar)dMbar oc M 6 ~*' 95 dM^ar at redshift z = 0, which is 
close to the slope n(M)dM oc M~ 2 dM for the number den- 
sity of collapsed haloes. This is the so-called "overcooling" 
problem (Cole 1991; Blanchard et al. 1992). 

The final radius of the cold gas in rotational equilib- 
rium is related to the initial radius by conservation of angu- 
lar momentum (Fall and Efstathiou 1980). Approximately, 
ro ~ Amin(rv,r cooi ), with the dimensionless spin param- 
eter A ee J|£| 1/2 G _1 M- 5/2 ~ 0.05 ± 0.03 (Barnes & Efs- 
tathiou 1987; Efstathiou et al. 1988; Zurek et al. 1988). Pre- 
vious studies only used the mean value of A in this formula. 
Hereafter, we introduce the A distribution from Barnes and 
Efstathiou (1987) model C0-4 (their fig. 11). According to a 
fit based on fig. 3 of Fall and Efstathiou (1980), the exponen- 
tial disk which forms from the dissipative collapse of the gas 
has a length scale ro ^ 1.26A 1 ' 17 min(rv, r coo i) and a radius 
including half the cold baryonic mass r 1 / 2 /ro = 1.68. A dy- 
namical time scale in the disk-like core is td yn = 2-!rri/ 2 /V c . 
It is important to note that only disks can form in this for- 
malism. The formation of elliptical galaxies (and of bulges of 
spiral galaxies) has to be explained by the merging of disks. 
Kauffmann et al. (1993) and Cole et al. (1994) showed that 
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this merging process can easily explain the current fraction 
of gE among bright galaxies. 

2.3 Star formation 

Locally, the SFR depends on numerous physical parame- 
ters. Nevertheless, phenomenological studies seem to show 
that, on galaxy scales, the SFR per unit surface density is 
proportional to the total gas surface density (neutral plus 
molecular) divided by the dynamical time scale of the disk 
(Kennicutt 1989, 1997). So we shall hereafter assume that 
the star formation time scale t+ is proportional to the dy- 
namical time scale of the disk td yn and we introduce a first 
efficiency factor (5. With = /3td yn , we take: 

M ga3 (t) 



SFR(t) = 



(2) 



Fig. 3 shows the predicted t* histogram compared to the 
histogram of "Roberts times" for a sample of 63 bright disk 
galaxies observed by Kennicutt et al. (1994). The Roberts 
time is defined as t R = (M H i + M H2 )/SFR(t ) where M H i 
is the gas mass in the HI phase measured from the 21 cm 
line, Mh 2 is the gas mass in the Hi phase measured from 
the CO line, and SFR(to) is the total star formation rate 
measured from the Ha line, under an assumption about the 
shape of the Initial Mass Function (hereafter this is Salpeter 
IMF). The interesting result is that the model correctly pre- 
dicts the shape and width of the histogram. We emphasize 
that this agreement is due to both the range of halo densities 
(scaling as (1 + z co u) 3 ) and dimensionless spin parameters 
A. Without the A scatter, the predicted distribution would 
be about three times too narrow. The average value of the 
observed histogram can be reproduced by taking /3 ~ 100 
for our SCDM model. 

Finally, for sake of simplicity, we use a Salpeter IMF 
with index x = 1.35. Stars have masses 0.1 < m < 120M Q . 
We also assume that the mass fraction blocked in dark ob- 
jects with masses below 0.1 Mq is negligible. 

2.4 Stellar feedback 

The explosion of massive stars can expel gas from the galax- 
ies and quench star formation, leading to a strong increase 
of the mass-to-luminosity ratios in small objects. Observa- 
tionally, HI holes and X-ray superbubbles are good evidence 
that such galactic winds are present in galaxies. The stellar 
feedback is introduced in a similar way by most authors, fol- 
lowing the original work by Dekel and Silk (1986). By equat- 
ing the gas binding energy to the thermal energy ejected by 
supernovae, one gets: 



. V 



M gas (t)(-^) 2 V C 2 = e / T t (t') V s N EsNdt', 



(3) 



where t/sn is the number of SNe per unit mass of stars, de- 
pending on the IMF. For our choice of IMF, n s n = 7.4 10~ 3 
Mg 1 . The output mechanical energy of a SN is Esn ~ 10 51 
erg. The escape velocity at radius r < rv in a singu- 
lar isothermal sphere truncated at radius rv is V eS c(r) = 
V2V C (1 -ln(r/rv)) 1/2 . The maximum of the A distribution 
corresponds to r-i/i/rv — 0.05, leading to V eS c/V c ~ 2.8. 
Since much of this energy is subsequently radiated away, 
we insert a second efficiency factor < e < 1. After some 



from log L IR -14 ji 
to log L„=6 




100 
A (in ,um) 



Figure 7. Model spectra in the IR and submm, for IR luminosi- 
ties 10 6 , 10 s , 10 10 , 10 12 , and 10 14 L bolQ . Emissivity index of big 
grains: m = 2 (dashes), m = 1.5 (solid lines), m = 1 (dotted 
lines) . 



algebra, the mass fraction of stars F* forming before the 
triggering of the galactic wind at time tw is given by: 



Fi, 



M*{t w ) 



M^tw) + M gas (t w ) 



(1 + {Vhot/Vc, 



(4) 



with a = 2 and V hot = (I4/14 sc )(2 ??S ]v£s]v) 1/2 e 1/2 = 
310 e 1/2 km s _1 for the Salpeter IMF. Nevertheless, there is 
much uncertainty on these parameters because of the cooling 
of supernova remnants before wind triggering. For instance, 
numerical simulations seem to suggest that e ~ 0.1 (Thorn- 
ton et al. 1997). Cole et al. (1994) introduced a fit based 
on SPH simulations of galaxy formation in which most of 
the feedback effect is due to momentum exchange rather 
than to ISM heating (Navarro and White 1993). For a typi- 
cal feedback parameter fv = 0.1, the numerical simulations 
can be fitted with a = 5 and Vhot = 130 km s _1 . We shall 
hereafter take these values as our standard parameters. The 
situation is still complicated by the existence of non-local 
feedback processes in addition to the local ones. Efstathiou 
(1992) and Blanchard et al. (1992) suggest that high-2 re- 
ionization of the IGM could prevent cooling in haloes with 
circular velocities below V equ ~ 20 to 50 km s^ 1 , and possi- 
bly as high as ~ (200) 1//3 V equ in case of adiabatic collapse. So 
it is very likely that the overall quenching of dwarf formation 
depends on redshift. Introducing this redshift dependence of 
dwarf formation partly alleviates the problem of the steep 
slope of the luminosity function (see e.g. Kauffmann et al. 
1993). So the situation appears to be very complicated, in 
the absence of a global theory of feedback processes. In the 
following, we shall simply model the combination of local 
and non-local feedbacks by introducing a simple (1 + z coU ) 
dependence for Vhot- 
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Figure 8. Model spectra with m = 1.5 superimposed to a com- 
pilation of data in the FIR and submm: Chini et al. 1986 (mean 
value), Stark et al. 1989, Ealcs et al. 1989, Carico et al. 1992, An- 
drcani and Franccschini 1992, Clements et al. 1993, Franceschini 
and Andreani 1995 (mean value). The dashes-and-dots show a 
typical spectrum of cold dust in the Milky Way (Reach et al. 
1995). Models and data are normalised at 100 fim. 



2.5 Spectral evolution of the stellar population 

A model of spectrophotometric evolution is used to compute 
the age dependence of the gas content, the UV to NIR spec- 
tra of the stellar populations, and the mass-to-luminosity 
ratios. The stars are placed on the Zero-Age Main Sequence 
of the HR diagram according to the IMF. The models use 
compilations of stellar evolutionary tracks taking into ac- 
count the various stages of stellar evolution in order to com- 
pute at each time step the distribution of stellar populations 
in the HR diagram. This distribution is combined with a 
library of stellar spectra and gives the synthetic spectrum 
F+\ . At the end of their lifetimes, stars die and return a frac- 
tion of their mass to the ISM. The model which is used here 
is described in Guiderdoni and Rocca-Volmerange (1987, 
1988), and Rocca-Volmerange and Guiderdoni (1988), and 
includes upgraded stellar tracks from Schaller et al. (1992) 
and Charbonnel et al. (1996). 

Panel a of fig. 4 shows the grid of models with various t* 
from 0.33 to 20 Gyr, which are introduced in order to follow 
the relative gas content of the galaxy g(t) = M gaB (t)/Mb ar . 
The heavy elements synthetised by stars are injected into 
the interstellar medium. The metallicity of the gas is esti- 
mated from the Instant Recycling Approximation Z g (t) = 
— j/zhi(?(t) with a yield yz depending on the choice of the 
IMF. Panel b shows the time variation of Z g (t) for the grid of 
models. The photometric properties can be computed after 
taking into account the intrinsic extinction (see the follow- 
ing section). As an example, panel c gives the face-on B — V 
colours. 



2.6 The slope of the luminosity function 

Before examining the modelling of the IR emission of galax- 
ies, we would like to quickly comment on the slopes of the 
_B-band luminosity and gas mass functions obtained by this 
simple version of the semi-analytic approach. 

In spite of the variety of SFR histories, and of the strong 
variation of the mass-to-luminosity ratios in time scales of 
a few Gyr, the shape of the luminosity function in the ab- 
sence of feedback processes is surprisingly similar to that of 
the baryonic mass function. For 4>{Lb)(ILb oc L b cILb, our 
model with j3 — 100 and the standard mass loss (a = 5 and 
Vhot = 130 km s _1 ) gives s — —1.4, whereas Loveday et 
al. (1992) find s — — 1. Such an uncomfortable situation is 
a robust result of all the recent attempts to model galaxy 
formation and evolution in a semi-analytic approach (see 
Kauffmann et al. 1994; Cole et al. 1994). Although it seems 
in disagreement with the nearby surveys (but other surveys 
seem to suggest an increase of the slope at the faint end, e.g. 
Marzke et al. 1994), this steep slope is necessary to repro- 
duce the deep redshift surveys (down to Bj < 24) and the 
faint galaxy counts (down to Bj < 28). Subtle selection ef- 
fects due to surface brightness could explain the discrepancy 
between the nearby luminosity function and the high-z one 
(McGaugh 1994; Lobo and Guiderdoni 1997). 

Indeed, this type of selection effect can be suspected be- 
cause the predicted gas mass function is in better agreement 
with the observational data. The predicted mass distribu- 
tion of gas at redshift 2 = has a slope n(M ga s)dM gas oc 
M g ~a a 3 dM gas (with the standard values or the parameters). 
Smaller objects form earlier on an average, with higher den- 
sities, smaller td yn , smaller t*, and, consequently, they are 
more affected by mass loss and their relative gas content 
is lower. It is possible to correct statistically the total gas 
mass function in order to predict an HI mass function and 
compare it to the observational slope —1.35 determined by 
Briggs and Rao (1993). In the sample of 63 disk galaxies 
gathered by Kennicutt et al. (1994), there is no systematic 
trend for Mhi /M gas versus M ga3 , ta, or the SFR. We have 
made an histogram of \og(MHi/M gas ), with values from —1 
to 0, and we have distributed the predicted number density 
at M g as into the range of Mhi values according to this his- 
togram. The effect of this correction turns out to be rather 
small. Thus there is a good agreement of the predictions 
with the observations. It is worthwhile to emphasize that 
the slope —1.35 of the data depends on the volume correc- 
tions of the survey for low-mass objects and is somewhat 
uncertain. 



3 SPECTRAL EVOLUTION OF DUST 
EMISSION 

3.1 Dust absorption 

Part of the energy released by stars is absorbed by dust 
and re-emitted in the IR and submm ranges. The deriva- 
tion of the IR/submm spectrum is a three-step process: i) 
computation of the optical thickness of the disks; ii) compu- 
tation of the amount of bolometric energy absorbed by dust; 
iii) computation of the spectral energy distribution of dust 
emission. The modelling of these steps is not an easy task 
since it requires addressing confused issues such as the chem- 
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Table 1. Scenarios of galaxy evolution 



Name 


fburst 


f quiescent 


% of ULIGs 


Line code 




(P = 10) 


09 = 100) 






Q 





1 


% 


dots and small dashes 


A 


0.04(1 + z coll f 


1 fburst 


% 


solid line 


B 


0.04(1 + z coll y> 


1 fburst 


5 % at all z co n 


dotted line 


C 


0.04(1 + z coll y> 


1 fburst 


90 % for z coH > 3.5 


long dashes 


D 


0.04(1 + z coll y> 


1 fburst 


15 % at all 2 co ;; 


short dashes 


E 


0.04(1 + z coll y 


1 fburst 


1-cxp -0.02(1 + z coU ) 2 


dots and long dashes 



ical evolution of the dust, and the geometrical distribution 
of dust relatively to stars. 

We assume that the gas is distributed in an exponential 
disk with truncation radius r g and mean H column density 
< N H (t) >= M gaa (t)/lAm H nr 2 g . The factor 1.4 accounts 
for the presence of helium. If r25 is the isophotal radius at 
25 mag arcsec~ 2 , the observations give r s /r25 — 1.6 (Bosma 
1981). The r-25 radii are consistently computed from the 
B magnitudes and ro radii of the galaxies. Fig. 5 shows 
the mean total and gas surface densities inside r g at red- 
shift 2 = 0. These surface densities fairly correspond to the 
crude estimate < Nfj(t) >= 6.8 W 21 g(t) atom cm~ 2 used 
in Guiderdoni and Rocca-Volmerange (1987) and Frances- 
chini et al. (1991, 1994). As noted by Guiderdoni & Rocca- 
Volmerange (1987), a galaxy with g ~ 0.20 has an H column 
density < Nh >— 1.4 10 21 atom cm~ 2 in good agreement 
with the observational value for late-type disks, in spite of 
the uncertainties in this estimate (Guiderdoni 1987). 

The mean optical thickness inside r g is given by: 



1.086 A v 

©V 7 I \ 



A v E B - 



,A X 



E B -v Nh 
< N H (t) > 
2.1 10 21 at cm- 2 



< N H (t) > 
)■ 



(5) 
(6) 



As in Guiderdoni and Rocca-Volmerange (1987) and 
Franceschini et al. (1991, 1994), the extinction curve de- 
pends on the gas metallicity Z g (t) according to power-law 
interpolations based on the Solar Neighbourhood and the 
Magellanic Clouds, with s = 1.35 for A < 2000 A and s = 1.6 
for A > 2000 A. The extinction curve for solar metallicity is 
taken from Mathis et al. (1983). 

It is not clear whether the disks of "normal" spirals 
are optically thin or optically thick. Catalogues of galaxies 
(such as the classical RC2 and the RSA) suggest t b = 0.7 
for spirals. On the other hand, several studies have pointed 
out the possibility of optically-thick disks (see e.g. Disney 
et al. 1989 and references therein). Since the gas decrease 
is counterbalanced by the metallicity increase, the optical 
thickness of our grid of models is maximal for g(t) = e~ s ~ 
0.20 with s = 1.6. Late-type disks spend most of their life 
time at optical depths tb — 0.4 - 0.7, in good agreement 
with what is suggested in the RC2 or RSA, and with the 
estimates of Rowan-Robinson (1992) based on IR A S results. 
Fig. 6 gives the predicted distribution of optical thicknesses 
which fairly compares with a sample of normal spirals and 
luminous IR galaxies (Franceschini and Andreani 1995). 

As in Guiderdoni and Rocca-Volmerange (1987) and 
Franceschini et al. (1991, 1994), wc also assume a simple 
geometric distribution where the gas and the stars which 



contribute mainly to dust heating are distributed with equal 
height scales in the disks. If t\ (t) is the optical thickness of 
the disks at wavelength A and time t, the extinction correc- 
tion (averaged over inclination angle i) is: 

A x {t) = -2.5 log < l-^(-™(t)/<*»i) (?) 
a\T\(t)/ cos i 

providing that the stars and gas have the same height scales. 



The factor oa 



(1 



L0\ 



,V2 



crudely takes into account 



the effect of the albedo u>\ (from Draine and Lee 1984) for 
isotropic scattering (Natta and Panagia 1984) and relates 
extinction to absorption. This "slab" geometry is interme- 
diate between the "screen" geometry exp — a\T\/ cos i and 
the "sandwich" geometry (with zero height scale for dust) 
(1+cxp — a\T\/ cos i) /2 which respectively lead to larger and 
smaller absorption. Observational analyses seem to suggest 
that this "slab" geometry is the best guess to fit the data (see 
e.g. Franceschini and Andreani 1995; Andreani and Frances- 
chini 1996). 

Finally, the bolometric luminosity which is absorbed by 
dust and released in the IR/submm is: 



Li R (t) = J F, A (t)(l-10- a4 ^ (i) )dA. 



(8) 



Panel c and d of fig. 4 respectively show the face-on extin- 
guished B — V colours and the IR bolometric luminosities 
per unit mass of galaxy for a grid of models with various 
t*. The corresponding range of IR-to-blue luminosity ra- 
tios 0.06 < Lib/ XbLb < 4 is very similar to the observed 
range for optically-selected samples, for instance the sample 
drawn from the Zwicky Catalogue with ra z < 15.7 mag in 
Soifer et al. (1987), and characterizes "mild starbursts" and 
LIGs. The IR-to-blue luminosity ratios are sensitive to tb 
and the geometry. If there are tidally-induced gas inflows 
resulting in smaller r g and to a "screen" -like gas distribu- 
tion, this ratio can easily reach 10-100 as in the ULIGs. We 
hereafter keep the LIG-type efficiency of conversion between 
the optical and the IR as a conservative estimate, and we 
shall consider the population of ULIGs in Sec. 4.3. 



3.2 Dust emission 

The emission spectra of galaxies are computed as a sum 
of various components, according to the method devel- 
opped by Maffei (1994), which uses the observational cor- 
relations of the IRAS flux ratios 12/im/60/im, 25pim/60^im 
and 100/im/60/im with Lir (Soifer and Neugebauer 1991). 
These correlations are extended to low Lir with the samples 
of Smith et al. (1987) and especially Rice et al. (1988). 
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Figure 9. Evolution of the cosmic comoving star formation rate 
density as computed from rest— frame UV luminosity densities, by 
using Salpctcr IMF with slope 1.35. The UV densities are based 
on the integration of the luminosity function fit on all magnitudes. 
The UV fluxes have not been corrected for intrinsic extinction. 
The solid squares are derived from the rest-frame 2800 A lumi- 
nosity density in the Canada-France Redshift Survey, and the 
local value is from the B-band luminosity function and average 
UV-.B colour (Lilly et al. 1996). The open triangles take into 
account NIR data to compute photometric redshifts in the Hub- 
ble Deep Field (Connolly et al. 1997). The solid triangles come 
from another photometric-redshift analysis of the HDF (Sawicki 
et al. 1997), while the open squares with the arrows are lower and 
upper values from the HDF with redshifts derived from Lyman— 
continuum drop-outs (Madau et al. 1996). Scenario Q ("quies- 
cent" mode, fi = 100) is plotted with dots and small dashes. It 
is unable to reproduce the steep decline of psfr between red- 
shifts 1 and 0. Other scenarios involve an increasing fraction of 
the "burst" mode (/3 = 10) with redshift, producing a population 
of LIGs. Scenario A (solid line) has no ULIGs. Various quanti- 
ties of ULIGs are included in scenarios B (dotted line), C (long 
dashes), D (short dashes) and E (dots and long dashes). See tab. 
1 for details and summary of lines codes. The SFRs slightly differ 
because t* is considered as an exponential time scale for the LIGs 
and a burst duration for the ULIGs. 



Several components are considered in the model spec- 
tra: 

• Polycyclic aromatic hydrocarbons (PAH). Because of 
their small size (< 1 nm), these molecules never reach ther- 
mal equilibrium when they are excited in a UV/ visible ra- 
diation field. Their temperature fluctuates and can reach a 
value much higher than the equilibrium temperature, ex- 
plaining the 12 /xm excess and the bands at 3.3, 6.2, 7.7, 
8.6 and 11.3 /jm. Their template emission spectrum is taken 
from the model by Desert et al. (1990). 

• Very small grains (VSG). They are made of graphite 
and silicates. These dust grains have sizes between 1 and 
10 nm. As the PAH, they never reach thermal equilibrium. 



Consequently, their emission spectrum is much broader than 
a modified black body spectrum at a single equilibrium tem- 
perature. The template emission spectrum is also taken from 
the model by Desert et al. (1990). 

• Big grains (BG). They are also made of graphite and 
silicates. These dust grains have sizes between 10 nm and 
0.1 fim. They (almost) reach thermal equilibrium and can be 
reasonably described by a modified black body £„-B„(Tbg) 
and emissivity e„ oc v m with 1 < m < 2. 

• Synchrotron radiation. This non-thermal emission is 
strongly correlated with stellar activity and, as a conse- 
quence, with IR luminosity (see e.g. Helou et al. 1985). 
According to observations at 1.4 GHz, this correlation is 
L„(1.4 GHz) = Li R /(3.75 10 12 x 10 9 ). Here L IR is in W, 
vso = 3.75 10 12 Hz is the frequency at 80 /u,m and q ~ 2.16 is 
determined from observations. Then we assume that we can 
extrapolate from 21 cm down to ~ 1 mm with a single aver- 
age slope 0.7, so that L v = L„(1.4 GHz)(y/lA GHz)' - 7 . 

The 60/100 colour gives the temperature Tbg of the 
BG, provided that an emissivity index m has been chosen. 
We hereafter test m — 1, 1.5 (standard value) and 2. Then 
the amount of BG, VSG and PAH are calculated iteratively 
from the 12/100, 25/100 and 60/100 ratios. The resulting 
spectra are computed from a few fim to several mm and 
evolve with Lm such as more luminous galaxies preferen- 
tially emit at shorter wavelengths. By construction, these 
spectra fit the IRAS colour correlations. Fig. 7 shows exam- 
ples of these model spectra for various Lib. and three values 
of the emissivity index for BG. 

In spite of its shortcomings, this method takes into 
account the observed spectral evolution in the IR and 
submm/mm ranges, while previous works generally use a 
single template spectrum in the whole spectral range. It 
is also possible to reproduce the FIR photometry of vari- 
ous individual galaxies by this method (Maflei 1994). Fig. 
8 shows a compilation of data superimposed to a series of 
model spectra. The compiled samples gather galaxies which 
have been observed at least at one wavelength in the submm 
range, (say, between 200 fim and 1.3 mm) with a positive de- 
tection. The various samples are very heterogeneous. They 
have been obtained with different telescopes and instru- 
ments and focus upon different types of galaxies ( "normal" 
spirals, LIGs, ULIGs). None of them is complete. Moreover, 
there is very little overlap between the samples, and the 
few galaxies which have been observed by various authors 
with various instruments generally have discrepant fluxes. 
These discrepancies can originate from the choices of beam 
size and beam shopping. Andreani and Franceschini (1992, 
1996) have studied the effect of the finite extension of the 
submm emission and computed an average aperture correc- 
tion, which turns out to be rather small. 

There is some debate about the presence of very cold 
dust which could emit at submm wavelengths. Chini et al. 
(1986), confirmed by Chini and Kriigel (1993), claim that 
the dust seen at IRAS wavelengths is unable to reproduce 
the observed submm emission. On the contrary, Stark et al. 
(1989), Eales et al. (1989), and Clements et al. (1993) do 
not detect any evidence for this new component. Chini et 
al. (1986) interpret their high fluxes as an evidence for the 
existence of cold dust. The value of the derived temperature 
depends on the choice of the emissivity index of these grains, 
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Figure 10. Evolution of the cold gas density parameter in 
damped Lyman-a absorbers. Solid squares: data without the 
APM QSO survey. Solid triangles: data including the APM QSO 
survey. Open squares: tentative correction for selection effects 
due to QSO obscuration (Storrie-Lombardi et al. f996). Solid 
hexagons: Natarajan and Pettini (1997). Open triangle: local es- 
timate from HI surveys (Briggs and Rao 1993). Line codes of 
scenarios Q, A, B, C, D, E are given in tab. 1. The various sce- 
narios involving the "burst" mode consume more gas than the 
"quiescent" mode (dots and short dashes). 

which can be anything between 1 and 2. Clearly, the uncer- 
tainties in the submm emission are large, and the need for a 
systematic survey of a complete sample is strong. There is 
no doubt that this will be one of the first targets of SCUBA. 
Because of these uncertainties, the models with an average 
index in — 1.5 plotted in fig. 8 seem to cover the range of 
observations at submm wavelengths, except the high fluxes 
observed by Chini et al. (1986). 



4 THE HISTORY OF STAR FORMATION IN 
THE UNIVERSE 

4.1 The "quiescent" mode of star formation 

In this section, we will introduce scenarios of evolution which 
will be used to compute the IR/submm properties of galax- 
ies, and to predict faint galaxy counts, and the CIB. The 
description of these scenarios and their line codes in the fig- 
ures are summarized in tab. 1. While a complete assessment 
of the energy budget of galaxies would require the moni- 
toring of multi-wavelength luminosity functions and galaxy 
counts (delayed to forthcoming studies), we hereafter only 
wish to address the issue of the overall evolution of the 
comoving SFR and gas densities in the universe. Pei and 
Fall (1995) and Fall et al. (1996) emphasized the relevance 
of these quantities to describe the overall evolution of the 
galaxy population. 



Figure 11. Rest-frame comoving luminosity density. Letters UV, 
B, I and IR respectively stand for 2800A, 4400 A, 10000 A and 
60 ^m. The emissivity at 1600 A is about 30 % higher than at 
2800 A. The luminosity densities arc based on the integration of 
the luminosity function fit on all magnitudes. Solid squares: lo- 
cal and Canada-France Rcdshift Survey (Lilly et al. 1996). Open 
triangles: NIR data are taken into account to compute photo- 
metric redshifts in the Hubble Deep Field (Connolly et al. 1997). 
Solid triangles: other estimates of photometric redshifts in the 
HDF (Sawicki et al. 1997). Open squares with arrows: HDF with 
redshifts from Lyman— continuum drop— outs (Madau et al. 1996). 
Solid hexagon: 60 fim local density corresponding to one third of 
the bolometric light radiated in the IR (Saunders et al. 1990). 
Line codes of scenarios A, B, C, D, E are given in tab. 1. The dif- 
ferent UV and IR emissions mainly result from different fractions 
of ULIGs (with a top— heavy IMF and strong extinction), with 
almost similar SFR histories. Strongly diffcrring high-2; IR emis- 
sion are obtained without being much constrained by the current 
status of (discrepant) UV/optical observations. 



Fig. 9 and 10 respectively show the predicted SFR and 
gas evolutions for scenario Q with j3 = 100. It corresponds 
to the fit of SFR time scales in disks, that is, the so-called 
"Roberts times" peaking at 3 Gyr and ranging between 0.3 
and 30 Gyr (Kennicutt et al. 1994). The feedback parame- 
ters are a = 5 and Vhot = 40(1 + z co u) km s _1 . This some- 
what arbitrary (1 + z) dependence is intended to mimic the 
effect of global re-heating at high redshift and avoids a step- 
wise change of Vhot at z — 2 caused by a bimodal regime. 
This has almost no effect at low z but it changes the high- 
redshift behaviour of the SFR density. 

Clearly the SFR density in scenario Q does not decline 
fast enough between z = 1 and the local universe, and, cor- 
respondingly, gas is not consumed sufficiently. The impossi- 
bility to reproduce the peak of SFR density with star forma- 
tion histories similar to those in disks is not surprising. We 
know that their present SFR time scales (the Roberts times) 
are large and that their SFRs have not changed significantly 
during the last few billion years (Kennicutt et al. 1994). 
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Figure 12. Predictions of the diffuse backgrounds in the 
FIR/submm and in the optical compared to the current up- 
per limits and detections. The solid squares show the level of 
CO-B-B/DIRBE residuals from Hauscr (1995). The similar shapes 
of the residuals and dark sky suggest that the subtraction of fore- 
grounds has been incomplete and that the plotted values are up- 
per limits. The thick solid line gives the CIB spectrum from the 
re-analysis of the COBE/FIRAS residuals, initiated in Puget et 
al. (1996) and revisited in Guiderdoni et al. (1997). The solid 
hexagons show the Cosmic Optical Background (COB) obtained 
by summing up faint galaxy counts down to the Hubble Deep 
Field limit. Strictly speaking, this is only a lower limit of the ac- 
tual COB, but the shallowing of the U and B-band counts sug- 
gests near— convergence at least at those wavelengths (Williams et 
al. 1996). The short dashes and long dashes give the prediction 
for no— evolution integrated up to redshift Zf or = 8 in a cosmology 
with h = 0.5 and Qq = 1. The other curves are computed for the 
SCDM model with h = 0.5, fio = 1, 08 = 0.67, for scenarios Q, 
A, and E plotted with the line codes of tab. 1. Scenarios Q and 
A arc not sufficient to reproduce the CIB. This suggests the ex- 
istence of an additional population of ULIGs taken into account 
in scenario E. 



The prediction for the evolution of the comoving SFR 
density is very similar to the result given by Baugh et al. 
(1997), though their model involves a more accurate pro- 
cedure to compute the merging history of haloes (Cole et 
al. 1994). As in our model, they choose a SFR proportional 
to the ratio of the cold gas content to the dynamical time. 
But galaxies within haloes can merge after spiralling to the 
centre of the haloes induced by dynamical friction. During 
the merging, all the available gas is quickly converted into 
stars. In spite of this explicit implementation of a kind of 
"burst" mode, starbursting does not seem to affect signifi- 
cantly the SFR history, and the evolution they predict does 
not reproduce the steep decline of the SFR density between 
z = 1 and 0. Probably the modelling of dynamical friction 
alone underestimates the true magnitude of the interaction 
and merging processes between galaxies in merging haloes. 



4.2 The "burst" mode of star formation 

We now wish to introduce another mode of star formation 
with /3 = 10. The SFR time scales are now much shorter. 
As a result, the evolution of the luminosity function differs 
from the evolution in scenario Q involving the "quiescent" 
mode of star formation. As shown in fig. 2, the luminosity 
function of the "quiescent" mode is built up by the accumu- 
lation of light coming from all galaxies, and its evolution 
reproduces the evolution of the mass distribution of col- 
lapsed haloes. On the contrary, the luminosity function of 
the "burst" mode reproduces the evolution in the formation 
of haloes. Each generation of halo formation is accompa- 
nied by a strong starburst which acts as a transient beacon. 
As a consequence, the evolutionary trends of the luminosity 
functions in the two modes strongly differ. 

Since we know that the local universe is dominated by 
the "quiescent" mode, we now consider a mix of two broad 
types of populations, one with a "quiescent" star forma- 
tion rate, the other proceeding in bursts. For the "burst" 
mode, we take an involved mass fraction increasing with 
redshift f bu rst(z) = ./w s t(0)(l + Zcoii) 1 , as suggested by the 
increasing fraction of blue objects showing tidal and merger 
features at larger z (Abraham et al. 1996). Such a depen- 
dence is also consistent with theoretical considerations on 
the merger rates of galaxy pairs in merged haloes (Carlberg 
1990). Wc limit the scope of the present paper to this phe- 
nomenological modelling of fburst(z). This is clearly a point 
which should be refined in forthcoming studies. Noting that 
the frequency of galaxy pairs is oc (1 + z) s with S between 2 
and 6 at ±la (Zepf and Koo 1989; Burkey et al. 1994; Carl- 
berg et al. 1994), we choose here a high evolution rate 7 = 5. 
Then /burst (0) is set to 0.04 in order to fit the SFR density 
at low z, resulting in an "all-burst" behaviour at z coH > 0.8. 
This will be our scenario A. As it appears from fig. 9, this 
phenomenological description of the increasing importance 
of bursts reproduces the steep decline of the SFR density. 
The origin of this fast evolution has still to be elucidated by a 
more exhaustive modelling of all interaction processes in the 
semi-analytic codes which follow the merging history trees 
of haloes and galaxies. Fig. 9 also shows that the predicted 
cosmic SFR history has a high-z evolution closer to the one 
derived from photometric estimates of redshifts rather than 
that from Lyman-continuum drop-outs. For instance, one 
should note that at z = 4, the predicted SFR is ten times 
the lower limit derived from the Lyman-continuum drop- 
outs without extinction. As a consequence of the high SFR, 
the gas content strongly evolves between z ~ 2 and 0, as 
shown in fig. 10. 

We can now compute the corresponding IR/submm 
emission by taking the conservative estimate of the average 
optical thickness and "slab" geometry as in the "quiescent" 
mode. As a result, this population of "mild starbursts" and 
"luminous UV/IR galaxies" (LIGs) have IR-to-blue lumi- 
nosity ratios in the range 0.06 < Lir/\bLb < 4 which is 
characteristic of blue-band selected samples (Soifer et al. 
1987), and should be fitted to the Canada-France Redshift 
Survey (selected in the observer-frame Iab band, roughly 
corresponding to the B band at z ~ 1), and to high-z HDF 
galaxies. Galaxies at the peak of the SFR time-scale distri- 
bution {U = 0.3 Gyr) have L IR /M ~ 6.3L(, o!Q /M©. Never- 
theless, their colours can still be very blue during the burst 
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Figure 13. Predictions of the diffuse backgrounds in the 
FIR/submm (blow up), compared to the acceptable range of the 
CIB at ±1(7 per point (Guiderdoni et al. 1997). Scenarios Q, A, B, 
C, D, and E are plotted with line codes of tab. 1. The predictions 
for the COB are similar to those in fig. 12. 

(B — V = 0.1 at 0.5 Gyr). The evolution of the comoving 
luminosity density in various UV/ visible bands and at 60 
/xm is compared to observational estimates in fig. 11. The 
local energy budget and its evolution from z = to 1 seem 
to be fairly reproduced. The difference between the fits of 
the comoving SFR and luminosity densities (fig. 9 and 11) 
originates in the presence of dust absorption. The local UV 
flux is correctly fitted in fig. 11, while the local SFR in fig. 
9 is about twice the value deduced from the local UV flux 
under the assumption of no extinction. 

By integrating the contribution of sources along the 
line of sight, we can predict the corresponding diffuse back- 
grounds. While the results for both scenarios fairly repro- 
duce the COB, as shown in fig. 12, scenario A does not do 
much better than scenario Q to reproduce the CIB. Their 
predictions are clearly barely compatible with the accept- 
able ±lcr range recalled in fig. 13. The mean amplitude of 
the CIB is twice the prediction, despite our high choice of 
7. Consequently, we can strongly suspect the existence of a 
population of galaxies which are more heavily extinguished. 

4.3 A heavily— extinguished component 

The above-mentionned CIB computed with scenario A 
seems to be something like a conservative estimate of the 
minimum IR/submm background due to typical CFRS and 
HDF galaxies. We now wish to assess how much star forma- 
tion might be hidden by dust shrouds and introduce an addi- 
tional population of heavily-extinguished bursts, which are 
similar to "ultra-luminous IR galaxies" , or ULIGs (Sanders 
and Mirabel 1996; Clements et al. 1996a). We maximize 
their IR luminosity by assuming that all the energy avail- 



Figure 14. Predicted 60 fim luminosity functions at z = for 
the SCDM model with h = 0.5, Qq = 1 and ag = 0.67. Scenarios 
Q, A, B, C, D, and E are plotted with line codes of tab. 1. At 
2 = 0, scenarios A, B, C, D, E give almost similar results since 
the "burst" mode only involves 4 % of the mass. Scenario Q gives 
a slightly higher luminosity function since more gas is left to fuel 
star formation at z = 0. Squares: observational luminosity func- 
tion for the IRAS Bright Galaxy Sample (Soifer and Neugebauer 
1991). Thick long dashes: observational luminosity function from 
a compilation of various samples (Saunders et al. 1990). 



able from stellar nucleosynthesis (0.007:rMc 2 ) is radiated in 
a heavily-extinguished medium, yielding a total IR lumi- 
nosity L IR /M ~ 42.5(a-/0.40)(t*/lGj/r)- 1 L 6oi0 /M , pro- 
vided the lifetimes of stars are smaller than the duration t t 
of the burst. We take < x >= 0.40 for stars with masses 
larger than ~ 5Mq (Schaller et al. 1992), and get luminosi- 
ties which are a factor ~ 20 larger than those of our mild- 
starburst/LIG mode. As a consequence of this simple model, 
the starburst only radiates in the IR/submm, and only dark 
remnants are left when it stops, without the slow evolution 
of low-mass stars. We distribute this population of ULIGs 
in two ways: i) A constant mass fraction of 5 % (scenario 
B) or 15 % (scenario D) at all z coU , mimicking a scenario 
of continuous bulge formation as the end-product of inter- 
action and merging, ii) 90 % of all galaxies forming at high 
ZcoU > 3.5 undergo a heavily-extinguished burst, mimick- 
ing a strong episode of bulge formation. These scenarios are 
now able to fit both the COB and CIB. Scenarios B and D 
seem more appropriate to reproducing the CIB at 300 /im 
while scenario C, with high-z ULIGs, has a stronger con- 
tribution at larger wavelengths. Of course, these last three 
cases are only illustrative, and a combination of these so- 
lutions would also fit the CIB. For instance, we introduce 
an ad hoc scenario E with a fraction of ULIGs increasing 
as 1 — cxp— 0.02(1 + z coH ) 2 . Such a dependence can be ob- 
tained if the fraction of ULIGs depends on the mean surface 
density and optical thickness of disks which roughly scale as 
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(1 + z co u) 2 in our modelling. As shown by fig. 12 and 13, 
scenario E nicely reproduces the COB and CIB, and could 
be considered as our "best fit" . 

It is clear from fig. 11 that none of the optical data re- 
flects the large differences between these scenarios, although 
the fraction of light in the IR varies widely at high z. In 
scenario A, the IR/UV ratio decreases with increasing z be- 
cause of the decreasing metallicity of galaxies. In scenario 
B and D, this effect in cancelled because the ULIG bursts 
are assumed to be optically-thick, with a top-heavy IMF, 
and the IR/UV ratio at z — 4 is similar to that at z = 0. In 
scenario C and E, the IR/UV ratio strongly increases with 
z and is ~ 100 times higher than for model B at z = 4. 
One should note that at this redshift, scenario C and E 
are roughly consistent with the lower limits of the UV- 
luminosity density derived from Lyman-continuum drop- 
outs, but with ten times as much SFR as directly derived 
if extinction is not taken into account. In these scenarios, 
galaxy formation at high z is an almost completely-obscured 
process. 

We emphasize that the family of scenarios summarized 
in tab. 1 is introduced within the same SCDM model. The 
dissipative and non-dissipative collapses are unchanged, and 
the characteristic time scale for the conversion of gas into 
stars is always proportional to the local dynamical time, 
which is the most natural time scale. Given that, we only 
change the "fuzzy" astrophysics introduced with the effi- 
ciency parameter (3. It seems reasonable to assume that the 
efficiency of star formation is typically an order of magni- 
tude greater for interacting galaxies (the "burst" mode) than 
for isolated galaxies (the "quiescent" mode), resulting in a 
(3 parameter an order of magnitude lower. While this as- 
sumption has yet to be put on firmer grounds, its effect on 
the global SFR and gas densities is very strong. As a sec- 
ond step, changes in the IMF and extinction can modify the 
optical and IR/submm energy budget. 



5 PREDICTIONS IN THE FIR AND SUBMM 

5.1 The FIR luminosity function 

Fig. 14 gives predictions for the z = luminosity function 
at 60 /jm compared with the observational determinations 
drawn from the IRAS Bright Galaxy Sample (Soifer and 
Neugebauer 1991) and from a compilation of various IRAS 
samples (Saunders et al. 1990). The difference of the two ob- 
servational luminosity functions at the faint end illustrates 
the uncertainties. The agreement of the predictions with the 
data is fair. It is not surprising that scenario Q seems to give 
a better fit, since most galaxies of the BGS are spirals, mild 
starbursts and LIGs. All the scenarios involving the burst 
mode (A to E) give similar results at z — since the fraction 
of mass involved in starbursts is low (4 %). Burst scenarios 
with shorter SFR time scales in the past consumed more gas 
than scenario Q and less fuel is now left for star formation. 

It has been recalled in Sec. 2.6 that this kind of semi- 
analytic models predicts too many low-luminosity galaxies 
in the optical bands at z = 0, with respect to the observa- 
tional field luminosity functions (see the references quoted in 
Sec. 2.6). As a consequence, the optical luminosity functions 
are too steep and can be reconciled with the observations by 




S (Jy) 



Figure 15. Predictions for differential galaxy counts (normalised 
to Euclidean counts) at 60 (im. Data is shown for IRA S counts at 
60 fim. Solid triangles: Faint Source Survey (Lonsdale et al. 1990). 
Solid squares: QMW survey (Rowan— Robinson et al. 1991b). 
Solid hexagons: revisited counts in the Very Faint Source Survey 
(Bcrtin et al. 1997). Open squares: North Ecliptic Pole Region 
(Hacking and Houck 1987). These latter counts might be affected 
by the presence of a super-cluster. The counts by Grcgorich et 
al. (1995), which arc plotted with open triangles, are probably 
contaminated by cirrus. The no-evolution model (for a cosmol- 
ogy with h = 0.5, f2o = 1) is shown with short dashes and long 
dashes. Scenarios A, B, C, D, E are plotted with line codes of 
tab. 1. Scenario D has too many local ULIGs and should be re- 
jected. Scenario E with an increasing fraction of ULIGs almost 
reproduces the flat behaviour of the counts. 

invoking subtle selection effects based on surface brightness. 
Our model also predicts too many low-luminosity galaxies 
in the IR, while the discrepancy is somewhat lower than in 
the blue. The slope is 4>{LiR)dLin oc Lj^' 25 dLjn. The nor- 
malisation at the bright end can be improved by a slight 
decrease (a few tenth dex) of the baryonic density flbar, or 
a slight increase of the baryonic mass-to-luminosity ratios, 
through the normalisation of the "dark mass" in the IMF. 
We do not attempt to apply this kind of fine tuning. 

5.2 Faint galaxy counts 

We hereafter explore the faint counts predicted by using the 
scenarios proposed in tab. 1. Fig. 15 gives the predictions 
for the IRAS 60 /im differential counts normalised to the 
Euclidean slope. Observational counts from the QMW sur- 
vey (Rowan-Robinson et al. 1991b), the Faint Source Survey 
(Lonsdale et al. 1990), the Very Faint Source Survey (Bertin 
et al. 1997), and the North Ecliptic Pole Region (Hacking 
and Houck 1987) are superimposed to the predictions. The 
latter survey might be partly affected by a large-scale struc- 
ture at z = 0.088 (Ashby et al. 1996). As shown by Bertin et 
al. (1997), the counts by Gregorich et al. (1995) are contam- 
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Figure 16. Predictions for differential galaxy counts (normalised 
to Euclidean counts) at 450 and 850 fim for observations with 
SCUBA. In contrast with counts at wavelengths smaller than 100 
fim, the submm counts are very sensitive to the details of the 
high-z evolution, because of the shift of the 100 fim bump into 
the observing bands. Scenarios A, B, C and E are plotted with 
line codes of tab. 1. 



inated by cirrus. The theoretical curves are not renormalised 
to the bright end of the observed counts. The scenarios in- 
volving an increasing fraction of the "burst" mode predict 
more galaxies than the no-evolution curve. Scenario D with 
15 % of ULIGs is rejected by the data. All the other scenarios 
have a moderate local fraction of ULIGs and are in agree- 
ment with the faint counts. Scenario E with an increasing 
fraction of ULIGs gives an almost flat curve which is rem- 
iniscent of the observational trend. Nevertheless, the rise 
of the counts below 0.1 Jy, if it is real, is not reproduced. 
However, our scenarios predict a correct amount of fluctu- 
ation and there is not much space for stronger evolution. 
More specifically, the counts predict a 60 /jm background 
fluctuation per beam in the Very Faint Source Survey (after 
removal of > 4a to t = 120 mjy sources) at the level of 14.1 
mJy (A), 16.0 mjy (B), 14.3 mjy (C) and 17.0 mjy (E), 
while the measured 68 % quantile is 30. 1± 1.2 mjy (Bertin 
et al. 1997). With a 25 mjy rms instrumental noise and 6.5 
mjy rms cirrus fluctuations (Gautier et al. 1992), there is 
still space for a 15. 41 2 ,; 2 , mjy fluctuation due to sources, in 
good agreement with our estimates. Moreover, if the cirrus 
fluctuations are non-gaussian, the rms value strongly over- 
estimates the 68 % quantile and can tolerate 16. 81?, 3 m Jy 
for source fluctuation. 

Fig. 16 and 17 give the predictions for various wave- 
lengths from 15 /im to 1.4 mm which correspond to current 
and forthcoming instruments. The IRAS 60 (jm counts are 
recalled in one of the panels, as well as the ISO-HDF 15 
fim deep counts (Oliver et al. 1997). Only scenarios A and 
E are plotted. Clearly there are three regimes: (i) Nearby 



galaxies are in the Euclidean zone, giving a count slope 
N(> S) oc S~ 3 ^ 2 . The value of the bright-end normalisa- 
tion at submm wavelengths depends on the choice of the 
emissivity index m. Since the count rates are very low, it 
is likely that we shall have to wait for the all-sky survey of 
PLANCK in order to fix the counts at the bright end. (ii) At 
short wavelengths, the curvature effect and the positive "k- 
correction" produce the bend of the faint counts; (iii) In the 
submm/mm range, the negative "k-correction" produces a 
bump in the faint counts, which reflects the passage of the 
100 fim emission bump into the observing bands. 

It appears from fig. 17 that the ISO-HDF counts sug- 
gest evolution, but that the predictions are not strongly sen- 
sitive to the details of the evolutionary scenario. Guiderdoni 
et al. (1997) have emphasized the relative degeneracy of the 
predictions at wavelengths shorter than 100 (im, including 
the IRAS 60 fim and the ISO 15 fim counts, and the strong 
sensitivity of the submm counts to the details of galaxy evo- 
lution, which is already apparent for ISO observations at 
175 /jm. Tab. 2 and 3 gather the count predictions for the 
various scenarios of tab. 1, and two typical sensitivity lev- 
els: 100 mjy which will be reached for point sources by the 
PLANCK all-sky survey, and which is comparable to the 
0.22 Jy level of the IRAS 60 fim Faint Source Catalogue; 
and 10 mjy which can be reached by a deep pencil-beam 
survey with SCUBA and FIRST. The number counts at 100 
mjy and 10 mjy are very sensitive to evolution and can 
differ by an order of magnitude at 350 fim, and more than 
two orders of magnitude at 850 /im. This explains the dis- 
crepancy of the predictions in the literature so far published. 
Hivon et al. (1997) will come back to the detection of submm 
sources and of the fluctuations they induce. As a conclusion 
of these predictions, while the IRAS data do not yield tight 
constraints on the evolution of galaxies at z > 0.2, the forth- 
coming deep survey with ISO at 175 /j,m and with SCUBA 
in the atmospheric windows at 450 and 850 fim (for which 
the array is matched to the instrument optics) will quickly 
help discriminate between the various scenarios, in the ex- 
pectation of PLANCK and FIRST. 

5.3 Redshift distributions 

Fig. 18 shows the prediction for the redshift distribution 
corresponding to the spectroscopic follow-up of two IRAS 
samples: the North Ecliptic Pole Region observed by Ashby 
et al. (1996) with 60 /urn fluxes 60 < Seo < 150 mjy, and the 
sample of Clements et al. (1996a) for galaxies of the Faint 
Source Catalogue with 5*60 > 300 mjy which are 'IR loud" 
(IR-to-blue ratios larger than 10). The data of the NEPR 
are almost complete and can be directly compared to our 
predictions. Nevertheless, the first two bins are probably af- 
fected by a super-cluster at 2 = 0.088. The redshifts of the 
"IR-loud" FSC sample have been taken by the authors in or- 
der to extract a subsample of ULIGs and not for a complete 
redshift survey. The histogram given for the galaxies which 
are not ULIGs (144 sources) has been roughly rescaled to 
take into account the number of sources not quoted in their 
Appendix because they already have redshifts in the liter- 
ature. The shaded histogram shows their 91 ULIGs. Then 
the histogram is renormalised in order to peak at the same 
relative level as the NEPR. The sample has only a 60 % 
completeness. In general, the scenarios seem to predict more 
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Figure 17. Predictions for number counts at various wavelengths for operating and forthcoming observing facilities and experiments: 
ISOCAM and ISOPHOT on-board ISO, SCUBA, FIRST, and PLANCK Kigh Frequency Instrument. Some of the curves are redundant. 
The curves correspond to the wavelengths from top to bottom, except in the IRAS/ ISO Panel where it is from bottom to top. Note 
that the 15 fim curve only takes into account dust emission (without the stellar component), and the observer— frame fluxes should be 
considered as lower values beyond z ~ 1. The no— evolution prediction at 15 fim (for h = 0.5 and Qq = 1) is plotted with short dashes 
and long dashes. Open stars: IRAS at 60 fim (Lonsdale et al. 1990). Solid hexagon: IRAS at 12 fim (Rush et al. 1993). Solid squares: 
ISO-HDF counts at 15 fim (Oliver et al. 1997). Scenarios A and E are plotted with line codes as in tab. 1. 



galaxies at high z than in the NEPR sample. The "IR-loud" 
FSC sample plotted here suffers from problems of overall 
normalisation and incompleteness, but it signals the pres- 
ence of galaxies at z > 0.2. The difference between the red- 
shift distributions for scenarios E and A gives the fraction 
of ULIGs in scenario E (there are no ULIGs in scenario A). 
The distribution of ULIGs with z seems to be much flatter 



than the observed peak in the "IR-loud" FSC sample. How- 
ever, it is difficult to assess the level of discrepancy with a 
60 % the incompleteness. Clearly the redshift follow-up of 
deep IR-selected samples is still an on-going project, and it 
will bring crucial information on the nature of the sources. 

Fig. 19 shows predictions for the redshift distributions 
at 60, 200, 350, and 550 fim, corresponding to various flux 
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Table 2. Predictions of galaxy counts in PLANCK HFI bands at 
the 100 mjy level (log number of sources sr _1 ). 



40 



Name 


350 (im 


550 (im 


850 fim 


1380 fim 


Q 


3.97 


3.06 


2.15 


1.18 


A 


3.93 


2.73 


1.67 


0.71 


B 


4.41 


3.48 


2.23 


0.83 


C 


4.95 


4.88 


4.39 


3.15 


D 


5.04 


4.34 


3.33 


1.77 


E 


5.08 


4.57 


3.83 


2.39 



Table 3. Predictions of galaxy counts in SCUBA and FIRST 
bands at the 10 mjy level (log number of sources sr _1 ). 



Name 


350 (im 


450 (im 


650 fim 


750 ^m 


850 (im 


Q 


5.66 


5.23 


4.41 


4.09 


3.80 


A 


5.98 


5.63 


4.73 


4.40 


3.98 


B 


6.27 


5.99 


5.34 


5.07 


4.75 


C 


6.80 


6.81 


6.64 


6.52 


6.37 


D 


6.66 


6.47 


6.01 


5.80 


5.55 


E 


6.84 


6.75 


6.45 


6.29 


6.09 



T 



data: 

solid lino: 60 SS 60 <150 mjy (NEPR) 
dashes: S e0 §300 mjy and IR loud 



upper' and lower curves: 
60 SS 60 S300 mJy 
300 <S 60 <IS00 mJy 




depths. Scenarios A and E (as well as B, C, and D not 
plotted here) give different predictions at z > 1. For scenario 
A, the redshift distributions at the 10 mjy flux level peak at 
z ~ 1, with a high-z tail encompassing 90 % of the sources 
at z ~ 2.5 to 3. Scenario E has galaxies at still higher z. In 
contrast with the IRAS samples at 60 (jm which, at the 100 
mjy flux level, only probe the universe at z ~ 0.2, the future 
spectroscopic follow-up of submm observations at the 10 
mjy level should be able to explore the IR/submm evolution 
of galaxies at z > 1. Finally, fig. 20 shows the contribution 
of sources fainter than S to the background value at 60, 200 
and 550 /im. The 10 mjy level which will be reachable by 
the forthcoming submm observations will allow the surveys 
to begin "breaking" the CIB into discrete units. 



6 DISCUSSION AND CONCLUSIONS 

We have used a simple semi-analytic model of galaxy forma- 
tion to derive the IR/submm statistical properties of galax- 
ies: luminosity function, faint galaxy counts, redshift dis- 
tributions, and the diffuse background. The model shares 
many common features with previous semi-analytic works 
focussing upon the optical properties of galaxies: (i) the 
collapse of the perturbations is described by the classical 
top-hat model under the assumptions of homogeneity and 
sphericity; in this simple version, we used the peaks formal- 
ism to compute the formation rate of haloes, (ii) the dis- 
sipative cooling and collapse is introduced, with the usual 
"overcooling" problem which can be partly suppressed by in- 
troducing feedback due to overall re-ionization and galactic 
winds triggered by SNe. (iii) star formation is proportional 
to the ratio of the gas content to the dynamical time scale 
of the galaxies; and (iv) stellar evolution is explicitly imple- 
mented. 

In order to make specific predictions for the IR/submm 
wavelength range, this model includes the following assump- 
tions: 

(i) An average optical depth of the disks is implemented 



Figure 18. Predicted redshift distribution of the North Ecliptic 
Pole Region observed by Ashby et al. 1996. The observed galaxies 
have 60 fim fluxes in the range 60 < Sao < 150 mjy. The first two 
bins might be affected by a super-cluster. The sample of Clements 
et al. (1996a) is also plotted for sake of illustration. The reader is 
referred to this paper for the description of the selection criteria. 
Roughly, the galaxies have Sao > 300 mjy and are "IR loud" 
(IR-to-blue ratio larger than 10). See text for the details. The 
shaded histogram show the ULIGs of this sample. Scenarios A 
and E arc plotted with line codes of tab. 1. Their normalisation 
is relative. 

and scales as the relative gas content and gas metallicity 
(estimated through the assumption of Instant Recycling Ap- 
proximation) ; no evolution of the dust composition is con- 
sidered; the variation of the extinction curve with metallic- 
ity is based on an interpolation of the Milky Way and the 
Magellanic Clouds which gives an extinction proportional to 

(ii) A simple "slab" geometry is assumed where dust and 
stars are mixed with equal height scales; 

(iii) Finally, the dust emission spectrum scales with the 
bolometric IR luminosity, following the observational rela- 
tionship of IRAS colours with luminosities. Such a depen- 
dence is obtained by changing the abundance of three com- 
ponents: PAH, very small grains and big grains (and the 
temperature of the big grains). 

We introduced two modes of star formation. In the "qui- 
escent" mode", the distribution of the characteristic time 
scales for star formation peaks at 3 Gyr, and ranges from 
0.3 to 30 Gyr. In the "burst mode", the distribution of t* 
peaks at 0.30 Gyr, and ranges from 0.03 to 3 Gyr. The "qui- 
escent" mode is unable to reproduce the rapid evolution of 
the cosmic SFR and gas densities. The steep decline of the 
SFR and gas content since z = 2 can be accommodated 
by assuming an increasing fraction of mass involved in the 
"burst" mode. Since one knows that the fraction of pecu- 
liar objects with signs of interaction/merging increases with 
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Figure 19. Predictions for redshift distributions at four wavelengths, and in different flux ranges. The flux bin widths are A log 5 = 1, 
and the central flux values (log 5 in Jy) are shown in each figure (from top to bottom). Scenarios A and E are plotted with line codes of 
tab. 1. 



the depth of the survey, we suspect that the bursting be- 
haviour of the CFRS and HDF galaxies is due to interactions 
and we introduce a phenomenological z dependence fitting 
the evolution of the pair rates. Since LIGs and ULIGs are 
respectively interacting and merging systems, it is reason- 
able to estimate that, reciprocally, CFRS and HDF galaxies 
with peculiar morphologies should emit in the IR. For these 
"mild starbursts" and "luminous IR/UV galaxies", we as- 
sume a conservative range of IR-to-blue luminosity ratios 
0.06 < Lir/\bLb < 4 computed from our standard as- 
sumptions on average optical thickness and geometry. Fi- 
nally, in order to reproduce the bright end of IR galaxies 
with Lir/\bLb > 10, we also introduced a population of 



heavily-extinguished galaxies with massive star formation, 
similar to ULIGs. 

With these ingredients, we designed a family of evolu- 
tionary scenarios which fit the local overall energy budget 
and its evolution to z ~ 1. These scenarios also fit the COB 
estimate computed from faint galaxy counts, and are con- 
sistent with the observational range for the CIB. So we feel 
confident that we roughly reproduce the local energy bud- 
get as well as a plausible range for the integrated budget 
along the line of sight. The scenarios predict UV and IR 
emissions strongly differring at high z. With the same SFR 
history, different UV fluxes can be obtained with different 
extinctions and IMF. In some of the scenarios, the UV fluxes 
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Figure 20. Contribution of sources fainter than S to the back- 
ground value at 60, 200 and 550 fim. Scenarios A and E are plotted 
with line codes of tab. 1. 

at z = 4 fit the current measurements (or lower values) of 
Madau et al. (1996), but with ten times more SFR than de- 
duced under the assumption of no extinction. Such a factor 
10, which includes the contribution of heavily-extinguished 
ULIGs, is much larger than the factor 3 corresponding to the 
correction for extinction derived by Pettini et al. (1997) for 
the HDF galaxies with Lyman-continuum drop-outs. This 
quantifies the strong warning that the deduction of the SFR 
density from UV fluxes misses the optically-dark side of 
galaxy formation. 

This family of scenarios was obtained within the context 
of the same cosmological model, with similar prescriptions 
for dissipative and non-dissipative collapses, star formation 
and feedback. The changes only affect the efficiency of star 
formation in a dynamical time scale (the /3 parameter), the 
IMF (Salpeter or only massive stars), and the extinction (av- 
erage extinction with "slab" geometry, or strong, "screen" - 
like extinction). We saw the extreme sensitivity of the pre- 
dictions to the parameter choice for this poorly constrained 
astrophysics. 

The scenarios are consistent with IRAS data. They give 
similar predictions for the IRAS 60 /im luminosity function 
and faint counts, provided the local fraction of ULIGs is 
lower than a few %. They are in rough agreement with the 
redshift surveys of IRAS samples, in spite of some uncer- 
tainty in the data. By construction, the IR spectra also fit 
the correlation of IR colours with total IR luminosity in 
the Bright Galaxy Sample. The selected scenarios are con- 
sistent with the ±lcr range for the observed CIB and they 
are used to "disentangle" the background into faint galaxy 
counts. This CIB is the first "post-IRAS' observation which 
helps constraining the high-z evolution of galaxies in the 
IR/submm and the energy budget integrated along the line- 
of-sight, before ISO results. In contrast with the counts at 



FIR wavelengths which are rather degenerate, the submm 
counts are very sensitive to the details of the evolutionary 
scenarios. Thus our study illustrates the importance of the 
submm range to constrain the evolution of the energy bud- 
get of galaxies. As emphasized by Guiderdoni et al. (1997), 
our understanding of the formation and evolution of galaxies 
is so far entirely based on surveys in the UV/ visible window. 
We are missing the part of star/galaxy formation hidden by 
dust. We need to open a second window in the IR/submm, 
which is being explored by ISO and SCUBA, and will be 
extensively studied by forthcoming telescopes and satellites 
such as SIRTF, SOFIA, FIRST and PLANCK. 

There are clearly many approximations and shortcom- 
ings in our approach, which could be listed as a research 
programme for forthcoming papers, (i) We need to imple- 
ment the astrophysical processes in a code which explicitly 
follows the merging history trees of haloes and galaxies. This 
code would have to include the various aspects of interac- 
tion between galaxies, and to explain why the interaction 
rate strongly declined since z = 1. It would also have to 
include the effect of interactions on the SFR. (ii) The mod- 
elling of the IR emission could take into account the evolu- 
tion of dust properties with chemical abundances, as well as 
more sophisticated transfer, (iii) The global energy budget 
of galaxies has to be addressed by fitting the faint galaxy 
counts in the UV/visible and K, as well as the forthcoming 
NIR counts from ISO, in particular at 15 /im. First results 
are already given here, and suggest an amount of evolution 
which is consistent with our predictions, (iv) The strong 
evolution of the comoving SFR density has consequences on 
the gas and heavy element abundances. The heavy elements 
made by the high SFRs are to be found in the stars, the cold 
gas, or the hot gas of the intergalactic medium. If the SFR 
turned out to be much higher than determined, for instance, 
by Madau et al. (1996), a significant fraction of the heavy 
elements should be present in the IGM. The amount of ele- 
ments synthetised in bursts would still be higher if the IMF 
is top-heavy, as in the simple model of ULIGs which is used 
here. According to Mushotzky and Loewenstein (1997), the 
absence of evolution of the iron content of clusters at z ~ 0.3 
is a first evidence that the global metal production is 2-5 
times higher than inferred from the UV drop-out technique. 

Despite the limitations listed above, this paper pro- 
vides the first predictions of the IR luminosity function, 
IR/submm galaxy counts, redshift distributions, and diffuse 
background, obtained from a semi-analytic model of galaxy 
formation and evolution which takes into account the main 
physical processes from the collapse of the density fluctua- 
tions to the absorption of UV and visible star light by dust 
and the re-emission in the IR/submm range. The predic- 
tions of this new model are consistent with the dynamical 
process of continuous galaxy formation predicted by the hi- 
erarchical growth of structures in a SCDM universe, and 
represent a significant progress with respect to previous phe- 
nomenological models designed to make predictions in the 
IR/submm. 

Tables giving the faint counts and redshift dis- 
tributions at various wavelengths between 15 /im and 
2 mm are available by anonymous ftp at ftp.iap.fr 
in the directory /pub/from_users/guider/fir, and on the 
Web node of the Institut d'Astrophysique de Paris at 
http: / / www.iap.fr / users / guider / fir .html . 
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APPENDIX A: MASS DISTRIBUTION OF 
COLLAPSED HALOES 

The initial, linear perturbation is assumed to be spherical 
and homogeneous. This is the so-called "top-hat" model, 
which has the friendly property of being entirely defined by 
two parameters, for instance the size R and density contrast 
(8p/p) z =o = So which are the linearly-extrapolated values at 
2 = 0, or, equivalently, by the mass M and collapse redshift 
Zcoii- If po is the current mass density of the universe, we 
have M = (4tt/3)R 3 p . 

The linearly-extrapolated density contrast, and the ra- 
tio of the radius of maximal expansion r m to the linearly- 
extrapolated size R, can be respectively computed as func- 
tions of the redshift of collapse: 



So = So[l + Z co ll], 



i/R 



i/-R[l + Zcol 



(Al) 



For instance, if Qo = 1, So = <5 C (1 + z co u) with 5 C = 1.68, 
and r m /R — (3/5So). After collapse and violent relaxation, 
a mean potential forms and the virialized halo ends as a 
singular, isothermal sphere truncated at "virial radius" rv — 
r m /2. If we define a "circular velocity" as V c = (GM/ry) 1 ^ 2 , 
the density profile of the relaxed halo at r < ry is: 



Ph{t) 



Vr 



(A2) 



and the mass included within radius r is simply M(< r) = 
Mr/ry- For instance, if Q.o — 1, the virial radius (in Mpc) 
and the circular velocity (in km s _1 ) are: 



rv = 0.17 ( 



V c = 160 ( 



M 



10 12 Af, 



o 



W3/i i \-l.-2/3 
") (1 + Zcoll) h 



Al 



10 12 M, 



a/3, 



1 



Zcoll 



) 1/2 h 1/3 . 







(A3) 



(A4) 



We hereafter take h = Ho/iWOkms' 1 Mpc' 1 ) 
virial temperature of the baryonic gas is: 



1 pm p 2 

2 k c 



35.9( 



Vc 



kms~ 



K, 



Finally, the 



(A5) 



m p is the mass of the proton, and the mean molecular weight 
/i = 0.59 corresponds to ionized gas with respective H and 
He mass fractions X — 0.75 and Y — 0.25. 

Then we compute the mass distribution of collapsed 
haloes from the peaks formalism (Bardeen et al. 1986), as 
in Lacey and Silk (1991) and Lacey et al. (1993). We start 
from the power spectrum of linear fluctuations P(k) which 
is normalized by as, the variance at a top-hat smoothing 
radius R — 8/i _1 Mpc. The peaks formalism uses mo- 
menta of the power spectrum with Gaussian smoothing 
on radius Rg- The mass included within such a radius is 
M = (27r) 3 / 2 R%po- So the relation between the Gaussian 
radius Rg and the top-hat radius R including mass M is 
simply Rg = 0.64_R. These momenta are computed from: 



<r m (R G ) 2 = g-j / k 2m P{k)\W k {kR G )\ 2 ^k 2 dk, 



(A6) 



where Wk(y) = exp(— y) is the Fourier transform of the 
Gaussian filtering function: 

1 l™|2 



W{x) 



(27^)3/2 



exp(- 



2R% 



(A7) 



The (comoving) number density distribution of peaks v 
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Sog I co(Rg) in which Sog is smoothed on the Gaussian scale 
Rg is: 



dv 
with: 



dv = n pk (R G )P(v, Rg)<1v, 



n pk (R G ) = 



1 V 

P(v,Rg) = — exp- — G(7,7i/), 



(A8) 

(A9) 
(AlO) 



such as J P(v, Ro)dv = 1. The functions R* and 7 are com- 
puted as: 



R*(R G ) = V3 



<ti{Rg) 



-t(Rg) 



<tI(Rg) 



(All) 



CT2(Rg)' (To(Rg)<J2{Rg) 

The constant c x — 0.6397 and an analytic fit of the function 
G(7, uj) is given by equ. (4.4) and (4.5) of Bardeen et al. 
(1986). If P(k) oc k n , it is easy to check that j(Rg) does 
not depend on Rg and that _R*(i?c) oc Rg- For £7o = 1, 
h = 0.5, and the CDM power spectrum with — 3 < n < 1, 
we have slow variations: 0.55 < ~((Rg) < 0.8 and 1.05 < 
R*(Rg) / Rg < 1.5 for 0.1 < R G < 100 Mpc. 

At that stage, we do not have the number density dis- 
tribution of peaks per bin of density contrast and radius. 
For that purpose, we have to introduce, as in Lacey & Silk 
(1991), a suggestion by Bond (1988) who interprets equ. 
(A9) as the total number density of peaks n pk (> Rg) at 
scale > Rg- By derivating equ. (A8), it is easy to get: 

d 2 n pk 



dRcdv 



dRadv - 



grfln^ 
Radln Rg 



n pk {R G )P{v, R G )dR G dv. (A12) 



We neglect the slow variation of 'y(Rg) with Rg- This is 
the number density of all peaks. We now need to count the 
peaks over a certain threshold v t h for collapse, and we get 
the formation rate of haloes: 



d 2 n for 
dRcdvth 



dRcdvth = 



d 
dv th 



roc 



d 2 n pk 
dRcdp 



dv)dR G dv t h- (A13) 



After a change of variables, we get the halo formation rate, 
since v t h = <5og,m>[1 + ZcoIi]/(To(Rg)- 



d 2 n fo 



d\nMd{l + z coU ) 



-n pfc (> R G ) 



din R* 



dlni?c 

dSoG.thU- + Zcoll 



' ct (Rg) d(l + Zcoii) 



P{y th ,Ra) 

(A14) 



In order to link the Gaussian smoothing to the top-hat for- 
malism, we take Rg = 0.64i? and SoG,th = Soao(R) / '(To(Rg) ■ 
The total number of collapsed peaks at redshift (1 + 2) is 
obtained by integrating equ. (A14) on all redshifts larger 
than (1 + z). 

The number density continuously increases with time, 
at all masses, and contrarily to what is found from the 
Press-Schechter formalism. The peaks formalism follows the 
collapse of all peaks, and counts high peaks as well as the 
broader, shallower peaks in which these high peaks are in- 
cluded. This so-called "cloud-in-cloud" problem results in 
a possible overestimate of the total number of galaxies. As 
a matter of fact, the total number density of peaks at the 
low-mass end of the mass function varies as M~ 2 , and the 
total mass density grows logarithmically with the lower mass 
cut-off. Anyhow, since other astrophysical phenomena will 



act to "suppress" galaxy formation in low-mass haloes, we 
find that this slow divergence is not too serious a problem. 

Finally, we do not explicitly include the merging of 
galaxies in merging haloes. In our crude modelling, new 
galaxies form at each generation of peak collapse. This model 
could be fitted to the description of "active" objects such 
as luminous galaxies seen in the IR. Such a simplifying 
assumption is very similar to that of Haehnelt and Rees 
(1993) who modelled QSO formation in each generation of 
halo formation. The issue of merging should be explicitly 
addressed by making Monte-Carlo realizations of the halo 
merging history tree and monitoring the merging of galax- 
ies in merged haloes, for instance with the dynamical fric- 
tion time scale (Kauffmann et al 1993; Cole et al. 1994; and 
following works). Such studies show that the galaxy merg- 
ing rate is relatively low, and that the resulting mass and 
luminosity functions do not differ significantly from those 
computed with cruder assumptions. More precisely, galaxy 
merging seems to be sufficient to make the 10 % fraction of 
giant galaxies which are elliptical, but does not produce a 
significant change in the slope of the luminosity function at 
faint magnitudes. The implementation of the astrophysics 
hereafter described into this type of code is clearly one of 
the next steps of our programme. Nevertheless, it is worth- 
while to note that only merging following the slow process of 
dynamical friction has been currently modelled. Two-body 
encounters and tidal interactions, which trigger starburst ac- 
tivity as it clearly appears from observational evidence, are 
not modelled yet in the current Monte-Carlo codes. 
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